#### prepare packages and functions ####
require(estimatr)
require(dichromat)

issue.var.name <- c("a.article9", "b.defense", "c.revisionism", "d.women", 
                    "e.gay", "f.immigrant", "g.growth", "h.tax")

color <- colorschemes$Categorical.12
plot.col <- c(color[12], color[12], color[12], 
              color[10], color[10], color[10], 
              color[2], color[2])

issue.label <- c("Revision of Article 9", "Increasing Defense Power", 
                 "Historical Revisionism", "Women's Empowerment", 
                 "Gay Marriage", "Accepting Foreign Workers", 
                 "Economic Growth", "Progressive Tax")

# function to compute MMs
MM <- function(data, attribute, outcome, id) {
  lm.data <- data.frame(id = data[, id], A = data[, attribute], 
                        Y = data[, outcome])
  l <- nlevels(lm.data$A)
  result.matrix <- matrix(NA, l, 3)
  for (i in 1:l) {
    level.label <- levels(lm.data$A)[i]
    result.matrix[i, ] <- unlist(lm_robust(Y ~ 1, data = lm.data, 
                                           subset = A == level.label, 
                                           clusters = id)[c(1, 6, 7)])
  }
  result.matrix
}

#### load data ####
task.data <- read.csv("task_data.csv")

# reorder the levels of attribute variables
position.label <- c("Agree", "Neither", "Disagree")
task.data$a.article9 <- factor(task.data$a.article9, levels = position.label)
task.data$b.defense <- factor(task.data$b.defense, levels = position.label)
task.data$c.revisionism <- factor(task.data$c.revisionism, levels = position.label)
task.data$d.women <- factor(task.data$d.women, levels = position.label)
task.data$e.gay <- factor(task.data$e.gay, levels = position.label)
task.data$f.immigrant <- factor(task.data$f.immigrant, levels = position.label)
task.data$g.growth <- factor(task.data$g.growth, levels = position.label)
task.data$h.tax <- factor(task.data$h.tax, levels = position.label)

#### analysis ####
## compute MMs
choice.result <- rating.result <- matrix(NA, 24, 3)
for (i in 1:8) {
  choice.result[(3 * (i - 1) + 1):(3 * i), ] <- MM(task.data, 
                                                   issue.var.name[i], 
                                                   "choice", 
                                                   "respondent.id")
  rating.result[(3 * (i - 1) + 1):(3 * i), ] <- MM(task.data, 
                                                   issue.var.name[i], 
                                                   "rating", 
                                                   "respondent.id")
}

# Figure 1
png("Figure_1.png", width = 5, height = 4.75, units = "in", pointsize = 9, 
    res = 2000, type = "cairo")
par(mar = c(3, 0, 0, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(4.58, 5.9), ylim = c(0, 33), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(seq(5.1, 5.9, 0.2), 0, seq(5.1, 5.9, 0.2), 33, 
         lwd = 0.4, col = "gray")
for (i in 1:8) {
  segments(5.08, 35 - 4 * i, 5.92, 35 - 4 * i, lty = 3, col = "gray")
  segments(5.08, 34 - 4 * i, 5.92, 34 - 4 * i, lty = 3, col = "gray")
  segments(5.08, 33 - 4 * i, 5.92, 33 - 4 * i, lty = 3, col = "gray")
  points(rating.result[(3 * (i - 1) + 1):(3 * i), 1], 
         c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i), 
         pch = 19, col = plot.col[i])
  segments(rating.result[(3 * (i - 1) + 1):(3 * i), 2], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i), 
           rating.result[(3 * (i - 1) + 1):(3 * i), 3], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i), col = plot.col[i])
  text(4.58, 36 - 4 * i, labels = issue.label[i], pos = 4)
  text(4.88, 35 - 4 * i, labels = "Agree", cex = 0.9, pos = 4)
  text(4.88, 34 - 4 * i, labels = "Neither", cex = 0.9, pos = 4)
  text(4.88, 33 - 4 * i, labels = "Disagree", cex = 0.9, pos = 4)
}
axis(1, at = seq(5.1, 5.9, 0.2), lwd = 0.5)
dev.off()

# Figure A.2
cairo_pdf("Figure_A2.pdf", width = 5, height = 4.75, pointsize = 9)
par(mar = c(3, 0, 0, 0), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.27, 0.6), ylim = c(0, 33), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(seq(0.4, 0.6, 0.05), 0, seq(0.4, 0.6, 0.05), 33, 
         lwd = 0.4, col = "gray")
for (i in 1:8) {
  segments(0.395, 35 - 4 * i, 0.605, 35 - 4 * i, lty = 3, col = "gray")
  segments(0.395, 34 - 4 * i, 0.605, 34 - 4 * i, lty = 3, col = "gray")
  segments(0.395, 33 - 4 * i, 0.605, 33 - 4 * i, lty = 3, col = "gray")
  points(choice.result[(3 * (i - 1) + 1):(3 * i), 1], 
         c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i), 
         pch = 19, col = plot.col[i])
  segments(choice.result[(3 * (i - 1) + 1):(3 * i), 2], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i), 
           choice.result[(3 * (i - 1) + 1):(3 * i), 3], 
           c(35 - 4 * i, 34 - 4 * i, 33 - 4 * i), col = plot.col[i])
  text(0.27, 36 - 4 * i, labels = issue.label[i], pos = 4)
  text(0.345, 35 - 4 * i, labels = "Agree", cex = 0.9, pos = 4)
  text(0.345, 34 - 4 * i, labels = "Neither", cex = 0.9, pos = 4)
  text(0.345, 33 - 4 * i, labels = "Disagree", cex = 0.9, pos = 4)
}
axis(1, at = seq(0.4, 0.6, 0.05), lwd = 0.5)
dev.off()

## comparison between choice-based and rating-based MMs
round(cor(choice.result[, 1], rating.result[, 1]), 3)

# Figure A.3
cairo_pdf("Figure_A3.pdf", width = 2.5, height = 2.5, pointsize = 8)
par(mar = c(4, 4, 1, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.4, 0.6), ylim = c(5.1, 5.9), 
     xlab = "Marginal Means for\nChoice-Based Outcomes", 
     ylab = "", xaxt = "n", yaxt = "n")
points(choice.result[, 1], rating.result[, 1], pch = 1)
axis(1, at = seq(0.4, 0.6, 0.05), lwd = 0.5)
axis(2, at = seq(5.1, 5.9, 0.2), lwd = 0.5)
mtext("Marginal Means for\nRating-Based Outcomes", side = 2, line = 2)
dev.off()

## test the significance of the differences in MMs between "agree" and "disagree" for economic issues
choice.AMCE.test <- lm_robust(choice ~ a.article9 + b.defense + 
                                c.revisionism + d.women + e.gay + 
                                f.immigrant + g.growth + h.tax, 
                              data = task.data, clusters = respondent.id)
summary(choice.AMCE.test)
rating.AMCE.test <- lm_robust(rating ~ a.article9 + b.defense + 
                                c.revisionism + d.women + e.gay + 
                                f.immigrant + g.growth + h.tax, 
                              data = task.data, clusters = respondent.id)
summary(rating.AMCE.test)